home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / zkscl.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  5.9 KB  |  153 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((zeror 0.0) (zeroi 0.0))
  12.   (declare (type double-float zeroi zeror))
  13.   (defun zkscl (zrr zri fnu n yr yi nz rzr rzi ascle tol elim)
  14.     (declare (type (simple-array double-float (*)) yr yi)
  15.              (type f2cl-lib:integer4 n nz)
  16.              (type double-float zrr zri fnu rzr rzi ascle tol elim))
  17.     (prog ((cyr (make-array 2 :element-type 'double-float))
  18.            (cyi (make-array 2 :element-type 'double-float)) (i 0) (ic 0)
  19.            (idum 0) (kk 0) (nn 0) (nw 0) (acs 0.0) (as 0.0) (cki 0.0) (ckr 0.0)
  20.            (csi 0.0) (csr 0.0) (fn 0.0) (str 0.0) (s1i 0.0) (s1r 0.0) (s2i 0.0)
  21.            (s2r 0.0) (zdr 0.0) (zdi 0.0) (celmr 0.0) (elm 0.0) (helim 0.0)
  22.            (alas 0.0))
  23.       (declare (type (simple-array double-float (2)) cyr cyi)
  24.                (type double-float alas helim elm celmr zdi zdr s2r s2i s1r s1i
  25.                 str fn csr csi ckr cki as acs)
  26.                (type f2cl-lib:integer4 nw nn kk idum ic i))
  27.       (setf nz 0)
  28.       (setf ic 0)
  29.       (setf nn (min (the f2cl-lib:integer4 2) (the f2cl-lib:integer4 n)))
  30.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  31.                     ((> i nn) nil)
  32.         (tagbody
  33.           (setf s1r (f2cl-lib:fref yr (i) ((1 n))))
  34.           (setf s1i (f2cl-lib:fref yi (i) ((1 n))))
  35.           (f2cl-lib:fset (f2cl-lib:fref cyr (i) ((1 2))) s1r)
  36.           (f2cl-lib:fset (f2cl-lib:fref cyi (i) ((1 2))) s1i)
  37.           (setf as (zabs s1r s1i))
  38.           (setf acs (- (f2cl-lib:flog as) zrr))
  39.           (setf nz (f2cl-lib:int-add nz 1))
  40.           (f2cl-lib:fset (f2cl-lib:fref yr (i) ((1 n))) zeror)
  41.           (f2cl-lib:fset (f2cl-lib:fref yi (i) ((1 n))) zeroi)
  42.           (if (< acs (- elim)) (go label10))
  43.           (multiple-value-bind
  44.               (var-0 var-1 var-2 var-3 var-4)
  45.               (zlog s1r s1i csr csi idum)
  46.             (declare (ignore var-0 var-1))
  47.             (setf csr var-2)
  48.             (setf csi var-3)
  49.             (setf idum var-4))
  50.           (setf csr (- csr zrr))
  51.           (setf csi (- csi zri))
  52.           (setf str (/ (exp csr) tol))
  53.           (setf csr (* str (cos csi)))
  54.           (setf csi (* str (sin csi)))
  55.           (multiple-value-bind
  56.               (var-0 var-1 var-2 var-3 var-4)
  57.               (zuchk csr csi nw ascle tol)
  58.             (declare (ignore var-0 var-1 var-3 var-4))
  59.             (setf nw var-2))
  60.           (if (/= nw 0) (go label10))
  61.           (f2cl-lib:fset (f2cl-lib:fref yr (i) ((1 n))) csr)
  62.           (f2cl-lib:fset (f2cl-lib:fref yi (i) ((1 n))) csi)
  63.           (setf ic i)
  64.           (setf nz (f2cl-lib:int-sub nz 1))
  65.          label10))
  66.       (if (= n 1) (go end_label))
  67.       (if (> ic 1) (go label20))
  68.       (f2cl-lib:fset (f2cl-lib:fref yr (1) ((1 n))) zeror)
  69.       (f2cl-lib:fset (f2cl-lib:fref yi (1) ((1 n))) zeroi)
  70.       (setf nz 2)
  71.      label20
  72.       (if (= n 2) (go end_label))
  73.       (if (= nz 0) (go end_label))
  74.       (setf fn (+ fnu 1.0))
  75.       (setf ckr (* fn rzr))
  76.       (setf cki (* fn rzi))
  77.       (setf s1r (f2cl-lib:fref cyr (1) ((1 2))))
  78.       (setf s1i (f2cl-lib:fref cyi (1) ((1 2))))
  79.       (setf s2r (f2cl-lib:fref cyr (2) ((1 2))))
  80.       (setf s2i (f2cl-lib:fref cyi (2) ((1 2))))
  81.       (setf helim (* 0.5 elim))
  82.       (setf elm (exp (- elim)))
  83.       (setf celmr elm)
  84.       (setf zdr zrr)
  85.       (setf zdi zri)
  86.       (f2cl-lib:fdo (i 3 (f2cl-lib:int-add i 1))
  87.                     ((> i n) nil)
  88.         (tagbody
  89.           (setf kk i)
  90.           (setf csr s2r)
  91.           (setf csi s2i)
  92.           (setf s2r (+ (- (* ckr csr) (* cki csi)) s1r))
  93.           (setf s2i (+ (* cki csr) (* ckr csi) s1i))
  94.           (setf s1r csr)
  95.           (setf s1i csi)
  96.           (setf ckr (+ ckr rzr))
  97.           (setf cki (+ cki rzi))
  98.           (setf as (zabs s2r s2i))
  99.           (setf alas (f2cl-lib:flog as))
  100.           (setf acs (- alas zdr))
  101.           (setf nz (f2cl-lib:int-add nz 1))
  102.           (f2cl-lib:fset (f2cl-lib:fref yr (i) ((1 n))) zeror)
  103.           (f2cl-lib:fset (f2cl-lib:fref yi (i) ((1 n))) zeroi)
  104.           (if (< acs (- elim)) (go label25))
  105.           (multiple-value-bind
  106.               (var-0 var-1 var-2 var-3 var-4)
  107.               (zlog s2r s2i csr csi idum)
  108.             (declare (ignore var-0 var-1))
  109.             (setf csr var-2)
  110.             (setf csi var-3)
  111.             (setf idum var-4))
  112.           (setf csr (- csr zdr))
  113.           (setf csi (- csi zdi))
  114.           (setf str (/ (exp csr) tol))
  115.           (setf csr (* str (cos csi)))
  116.           (setf csi (* str (sin csi)))
  117.           (multiple-value-bind
  118.               (var-0 var-1 var-2 var-3 var-4)
  119.               (zuchk csr csi nw ascle tol)
  120.             (declare (ignore var-0 var-1 var-3 var-4))
  121.             (setf nw var-2))
  122.           (if (/= nw 0) (go label25))
  123.           (f2cl-lib:fset (f2cl-lib:fref yr (i) ((1 n))) csr)
  124.           (f2cl-lib:fset (f2cl-lib:fref yi (i) ((1 n))) csi)
  125.           (setf nz (f2cl-lib:int-sub nz 1))
  126.           (if (= ic (f2cl-lib:int-sub kk 1)) (go label40))
  127.           (setf ic kk)
  128.           (go label30)
  129.          label25
  130.           (if (< alas helim) (go label30))
  131.           (setf zdr (- zdr elim))
  132.           (setf s1r (* s1r celmr))
  133.           (setf s1i (* s1i celmr))
  134.           (setf s2r (* s2r celmr))
  135.           (setf s2i (* s2i celmr))
  136.          label30))
  137.       (setf nz n)
  138.       (if (= ic n) (setf nz (f2cl-lib:int-sub n 1)))
  139.       (go label45)
  140.      label40
  141.       (setf nz (f2cl-lib:int-sub kk 2))
  142.      label45
  143.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  144.                     ((> i nz) nil)
  145.         (tagbody
  146.           (f2cl-lib:fset (f2cl-lib:fref yr (i) ((1 n))) zeror)
  147.           (f2cl-lib:fset (f2cl-lib:fref yi (i) ((1 n))) zeroi)
  148.          label50))
  149.       (go end_label)
  150.      end_label
  151.       (return (values nil nil nil nil nil nil nz nil nil nil nil nil)))))
  152.  
  153.